home *** CD-ROM | disk | FTP | other *** search
/ CD ROM Paradise Collection 4 / CD ROM Paradise Collection 4 1995 Nov.iso / science / earthqua.zip / EQSIM / GENDISP.CPP < prev    next >
C/C++ Source or Header  |  1995-02-04  |  4KB  |  138 lines
  1. // GENDISP.CPP - Generates the earthquake displacements to be used.
  2. // Created by Misha Koshelev.
  3.  
  4. #include <stdio.h>       // Standard functions
  5. #include <stdlib.h>      // Standard functions
  6. #include <math.h>                // Math stuff
  7. #include <time.h>                // For random numbers
  8. #include <alloc.h>               // For allocating memory
  9. #include <graphics.h>            // For error messages
  10. #include <dos.h>
  11. #include "gendisp.h"
  12.  
  13.                  // Global variables
  14. unsigned int waveletn;           // Number of wavelets
  15. double maxt;                     // Maximum parameter t in call to disp(t)
  16. unsigned int fract_dim;          // The fractal dimension
  17. double intensity;                // The Earthquake Intensity
  18. double far *offset = NULL;       // Array of random phases (or offsets)
  19. double far *steps = NULL;        // Steps
  20. double far *spectrand = NULL;    // Spectrum random numbers
  21.  
  22. // Inits global variables
  23. // Input: Intensity, Wavelet Number, Maximum time, the
  24. // step, and the fractal dimension.
  25. //
  26. // Maximum time is the biggest parameter time with which disp() is called.
  27. // Needed for a non-repetative wave.
  28. //
  29. // Output: None
  30. //
  31. void init_globals(double intens, unsigned int wavelet_n, double max_t,
  32.      double stp, double fd)
  33. {
  34.    if (stp >= 2 * M_PI / max_t)
  35.    {
  36.       closegraph();
  37.       printf("                          2 * PI\n");
  38.       printf("Step must be smaller than ------\n");
  39.       printf("                          %f    \n", max_t);
  40.       exit(1);
  41.    }
  42.    waveletn = wavelet_n;            // Set values
  43.    maxt = max_t;
  44.    intensity = intens;
  45.    fract_dim = fd;
  46.    if (offset && steps && spectrand)       // If arrays not empty destroy
  47.       done_globals();
  48.    offset = (double far *)farcalloc(waveletn, sizeof(double));
  49.    steps = (double far *)farcalloc(waveletn, sizeof(double));
  50.    spectrand = (double far *)farcalloc(waveletn, sizeof(double));
  51.    if (!offset || !steps || !spectrand)
  52.    {
  53.       closegraph();
  54.       printf("Not enough memory.\n");
  55.       exit(1);
  56.    }
  57.    for (int n=0; n<waveletn; n++)
  58.    {
  59.       offset[n] = dblrand(2 * M_PI, 4);
  60.       spectrand[n] = dblrand(1, 4);
  61.       steps[n] = stp * (n + 1);
  62.    }
  63. };
  64.  
  65. // Done globals
  66. void done_globals(void)
  67. {
  68.    farfree(offset);
  69.    farfree(steps);
  70.    farfree(spectrand);
  71.    offset = NULL; steps = NULL; spectrand = NULL;
  72. };
  73.  
  74. // A doubleing point random number generator
  75. // Input: Limit, Number of digits
  76. //
  77. // Output: Random number from 0 to limit
  78. //
  79. double dblrand(double limit, unsigned char digitn)
  80. {
  81.     int n;
  82.     time_t t;
  83.  
  84.     srand((unsigned) time(&t));                // Start random number generator
  85.     n = random((int)(limit * pow10(digitn)));  // Get a random integer
  86.     if (n < 0)                           // If that number is negative
  87.        n *= -1;                          // make it positive
  88.     return (n / pow10(digitn));   // Make integer a decimal number
  89. };
  90.  
  91. // Spectrum function.
  92. // Input: Omega, Earthquake Intensity
  93. //
  94. // Output:   I
  95. //         ----- * random
  96. //           fd
  97. //          w
  98. //
  99. double spectrum(double omega, double intensity)
  100. {
  101.     int idx = (int)(omega / steps[0]);
  102.  
  103.     return ((intensity / pow(omega, fract_dim)) * spectrand[idx]);
  104. }
  105.  
  106. // Generates the displacement for a certain time.
  107. // Input: Time
  108. //
  109. // Output: A displacement
  110. // Each element is:
  111. //        N
  112. //        --
  113. // d(t) = \   A(w ) sin (w t + l )
  114. //        /      n        n     n
  115. //        --
  116. //        n=0
  117. //
  118. // where d is the displacement, w  is n * step, the function A is the
  119. //                               n
  120. // spectrum function, and l is an array of random offsets.
  121. //
  122. double disp(double time)
  123. {
  124.     double cur_disp = 0.0;
  125.  
  126.     if (time > maxt)
  127.     {
  128.        closegraph();
  129.        printf("The parameter time as passed to function disp(time) exceeds\n");
  130.        printf("the maximum specified.\n");
  131.        exit(1);
  132.     }
  133.     for (int n=0; n<waveletn; n++)       // Loop wavelets
  134.      cur_disp += spectrum(steps[n], intensity) * sin((steps[n] * time) +
  135.             offset[n]);
  136.     return (double)cur_disp;
  137. }
  138.